Excitations and quantum fluctuations in site diluted two-dimensional antiferromagnets 
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We study the effect of site dilution and quantum fluctuations in an antiferromagnetic spin system 
on a square lattice within the linear spin- wave approximation. By performing numerical diagonaliza- 
tion in real space and finite-size scaling, we characterize the nature of the low-energy spin excitations 
for different dilution fractions up to the classical percolation threshold. We find nontrivial signa- 
tures of fractonlike excitations at high frequencies. Our simulations also confirm the existence of 
an upper bound for the amount of quantum fiuctuations in the ground state of the system, leading 
to the persistence of long-range order up to the percolation threshold. This result is in agreement 
with recent neutron-scattering experimental data and quantum Monte Carlo numerical calculations. 
We also show that the absence of a quantum critical point below the classical percolation threshold 
holds for a large class of systems whose Hamiltonians can be mapped onto a system of coupled 
noninteracting massless bosons. 

PACS numbers: 75.10.Jm, 75.50.Ee, 75.30.Ds, 75.40.Mg 



I. INTRODUCTION 



The problem of the interplay of quantum fluctua- 
tions and disorder in low dimensional systems is of 
fundamental importance in modern condensed matter 
physics. It is relevant for the understanding of the metal- 
insulator transition in metal-oxide-semiconductor field- 
effect transistors (MOSFETS)fi impurity effects in d- 
wave superconductors,^ non- Fermi-liquid behavior in U 
and Ce intermetallics,.? and the persistence of long-range 
order (LRO) in two-dimensional (2D) spin-1/2 quantum 
antiferromagnets (AFM))4 among others. 

The 2D square lattice with nearest-neighbor hopping 
undergoes a classical percolation transition upon ran- 
dom dilution. For the case of site dilution, the transi- 
tion occurs at the hole concentration Xc ~ 0.41^5. while 
for bond dilution the largest (infinite) connected clus- 
ter disappears when exactly half of the bonds are absent 
(i.e., Xc = 1/2).^'^ Thus, no ground-state long-range or- 
der is possible for any model with short range interac- 
tions in these lattices above Xc- For those models where 
order does exist in the clean limit, it is natural to ask 
whether dilution can enhance quantum fluctuations to 
the point of destroying long-range order at some dop- 
ing fraction below Xc- This possibility has led to several 
theoretical and experimental investigations in a variety 
of systems i^iSi^iiSiii In particular, a recent neutron scat- 
tering experimeniji in the site-diluted 5=1/2 Heisen- 
berg quantum AFM La2Cui_a:(Zn,Mg)j;04 (LCMO), in- 
dicates that LRO exists up to Xc- This fact was cor- 
roborated by an extensive quantum Monte Carlo (QMC) 
simulation* and by spin-wave theory (SWT) analytical 
calculation in the dilute regime^ The numerical and ex- 
perimental data suggest that the disappearance of or- 
der in the ground state is dominated by a classical effect 



and no quantum phase transition takes place below Xc- 
The analytical calculation,^ on the other hand, points to 
a nontrivial dilution-induced softening of low-frequency 
spin-wave excitations. However, the latter was carried 
out within the T-matrix approximation, which excludes 
coherent superposition and interference, and thus could 
not account for localization effects. 

In this work we carry out exact real space numerical 
diagonalizations of the dilute 2D Heisenberg AFM within 
the linear spin-wave approximation for dilution fractions 
ranging from the clean limit to the the classical percola- 
tion threshold. Although our method cannot be used to 
investigate systems as large as those used in previous nu- 
merical studies of the dynamical structure factor alon^ 
it provides complete access to eigenergies and eigenvec- 
tors, thus allowing us to probe more carefully into the 
structure of the excited states of the site dilute AFM. 
We find that excitations in this system break up into 
different modes as the amount of dilution is increased. 
The multipeak structure of the spectral function shows 
the simultaneous existence of extended (magnons) and 
localized (fractons) excitations, similarly to what was 
observed experimentally by Uemura and Birgeneau for 
dilute three-dimensional AFMs several years ago^ 

We also argue that the absence of a quantum critical 
point in diluted systems below universal feature 

for a large class of continuous models of the Heisenberg 
type that can be mapped into a system of coupled har- 
monic oscillators within some approximate scheme. We 
establish an upper bound for the amount of quantum 
fluctuations in these diluted noninteracting bosonic sys- 
tems and show that quantum fluctuations are bounded 
even at the percolation threshold for 2D systems. Indeed, 
one may wonder whether this generic behavior of diluted 
bosonic systems is related to universality classes dictated 
solely by symmetries of the disordered Hamiltonian, as 
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in the case of fermionic models. 

The paper is organized as follows. In Sec. ^Jwe de- 
fine the system Hamiltonian and find the corresponding 
non-Hermitian eigenvalue matrix problem in lattice co- 
ordinates within the linear spin-wave approximation. In 
Sec. ?? we present and discuss the numerical methods 
used to generate the random dilution, the identification 
of the largest connected cluster, and the diagonalization 
of the eigenvalue problem. The numerical results for the 
spin-wave excitation spectrum and the quantum fluctua- 
tion corrections to the AFM ground state are presented 
and discussed in Sec. ??. In Sec. ?? we argue that the 
absence of quantum phase transitions in site dilute sys- 
tems occurs for spin Hamiltonians that can be mapped 
onto a system of coupled harmonic oscillators. Finally, 
in Sec. ??, we draw our conclusions and point to future 
directions. 



II. DILUTE HEISENBERG 
ANTIFERROMAGNET IN THE SPIN WAVE 
APPROXIMATION 



We begin by reviewing the well-known connection be- 
tween magnetic and bosonic systems. In particular, we 
are interested in spin Hamiltonians of the form 



(1) 



where S"" is the a = x,y,z component of a spin S at 
site i, J a is the nearest-neighbor exchange constant, and 
rji = 0(1) if the site is empty (occupied). The empty 
sites are randomly distributed over the whole sample with 
uniform probability. Calling N the total number of sites, 
we define the fraction of occupied sites as 
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The dilution fraction is then defined as 

The isotropic AFM Heisenberg model corresponds to 
Ja = J > Q ioi a — x,y, z. When LRO is present 
(say, along the z direction), the spin operators can be 
written in terms of bosonic operators using the Holstein- 
Primakoff methodii^ Since the square lattice is bipartite, 
it can be divided up into two square sublattices, A = C^- 
and B — Thus, 
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with j G B. The bosonic operators obey the usual com- 
mutation relations, namely. 



[a^^a],] = 6u' and [bj,b],] 



(9) 



with all other commutators equal to zero. 

For the large spin case (5* 3> 1) or when the number of 
spin waves is small (n^ = (oiaj), rtj = {bj 6]) <C S), we 
can expand the square root in Eqs. (PJ and © in powers 
of Hi and rij, keeping only linear terms. That allows us 
to write the approximate bilinear bosonic Hamiltonian 
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The first term on the right-hand side of Eq. (|10|l repre- 
sents the classical ground-state energy of the AFM. Using 
the bosonic commutation relations, we can rearrange the 
Hamiltonian in the following form: 
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-JS{S- 



Hsw, 



(11) 



where, now, the first term on the right-hand side repre- 
sents the ground state energy in the absence of quantum 
fiuctuations, while the latter is described by the second 
term. 



Hsw = ^^^mTl3 {a\c 

(v) 
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Hereafter we will drop the constant ground state energy 
term and will only study the eigenstates of i?sw- The 
spin-wave Hamiltonian contains bilinear crossed terms 
which, separately, do not conserve particle number. At 
this point, it is worth simplifying the notation by drop- 
ping the distinction between bosonic operators living on 
different sublattices and reordering the summation over 
sites, 



(13) 



where both indexes in the sum run over all sites in the 
lattice. The matrices K and A are defined as 



(14) 
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(the sum run over all nearest-neighbor sites to i) and 
A / ^'^^ nearest neighbors 



0, otherwise. 



(15) 



Notice that both matrices K and A are real and sym- 
metric. 



A. Bogoliubov transformation 

It is possible to diagonalize the spin-wave Hamiltonian 
through an operator transformation of the Bogoliubov 
type, 



n 
n 

or, in matrix notation. 
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V* U* 
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where the column vectors a and a contain the operators 
Ui and an, respectively, while the N x N matrices U and 
V contain the coefRcients {ui„} and {um}, respectively, 
with i,n = 1, . . . , iV. Assuming that the new operators 
also obey the canonical commutation relations, 

[an,aln] = Sn7n and [an, am] = [a]^,aj„] = 0, (19) 

one arrives at the following constraints for the transfor- 
mation coefficients: 
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^ ^ i'^in '^jn '^in'Vjn) — ^i. 
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^ ^ (^m ^jn ^in ^jn) — 0- 
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In matrix notation. 



UU^ - VV'^ ^ I. 



N 
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where In is the N x N unit matrix. These relations can 
be put into a more compact form by defining the matrices 



and 



T = 



E = 



U V 

V* u* 
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Thus, Eqs. and become one, 

T E rt = E. (26) 
Since E^ — I2N, we find, after a simple algebra, that 

T^E T = E. (27) 

As a result we have two additional (though not indepen- 
dent) sets of orthogonality equations. 
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^{U*n Uim - Vin <„J = 5n 



and 
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^ ^ (^m Vim Vin Wj„) — 0. 
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B. Non-Hermitian eigenvalue problem (Ref. |l, 



The transformation defined by Eq. (|18|l allows us to 
diagonalize the spin-wave Hamiltonian in terms of the 
new bosonic operators. For that purpose, we choose T 
such that 



Q+ 





(30) 



where n± are diagonal matrices containing the eigenfre- 

quencies: [r2±]n„ = to^^ for n = 1,.. .,N. The eigen- 
value problem defined by Eq. (|30|l can be further simpli- 
fied. Recalling Eq. (|27|l . we have that 
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(31) 



In fact, it is not difhcult to prove that eigenfrequency 
matrices obey the relation fl^ = 171 = ^, with being 
a diagonal matrix with real entries only, as one physically 
expects. As a result, 

-V* u* I [ n 

(32) 

We can break up this 2N x 2N matrix equation into two 
coupled N X N matrix equations. 
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Ku + Av* = un, 
AU + KV* = -v*n, 



(33) 



or, alternatively, writing explicitly the matrix elements, 

^[KijUjn+ AijVjn] = (^nUin, (34) 
[Aij Ujn + K,^ V*^] = -U}nV*n, (35) 
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for all n and i. Thus, for a given eigenstate n, we can 
define an eigenvalue matrix equation in the usual form, 
namely. 



K A 
-A -K 



(36) 



(Notice that each u„ and w„ is now a column vector 
with components running through all i = 1, . . . , N lat- 
tice sites.) The 27V x 2N matrix shown in Eq. 
is clearly non-Hcrmitian, but its eigenvalues are all real. 
Notice also that if a;„ is an eigenvalue with corresponding 
eigenvector (u„,z;*), then —Un is also an eigenvalue, but 
with (w* , Un) as the corresponding eigenvector. Thus, de- 
spite the fact that the non-Hermitian matrix provides 2N 
eigenvalues (eigenfrequcncies) , we should only keep those 
N that are positive and whose corresponding coefficients 
M„ and Vin satisfy Eqs. ((201, (EU, I2H1), and if^ . 

The non-Hermitian matrix in Eq. H36|l contains only 
integer elements: 0, 1, 2, or 4 in the diagonal (correspond- 
ing to Kii, i.e., the number of nearest neighbors to site i) 
and or 1 in the off-diagonal components (correspond- 
ing to Aij, i.e., 1 when i and j are nearest neighbors and 
zero otherwise). It is strongly sparse, although without 
any particularly simple pattern due to the presence of 
dilution disorder. 

It is easy to verify that there are at least two zero 
modes in Eq. H36(l . i.e., two distinct nontrivial solutions 
with zero eigenvalue: 
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for all i 



. , N, and 



average staggered spin per site along the z direction can 
be written as 



TVfstagg 

niz = — — = S - Sm^, 



with 



Sm, = — 
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and 



(5mf = (a|ai). 
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Notice that brriz describes the spin-wave correction to the 
average staggered magnetization (always a reduction) . 

In order to express in terms of the coefficients 
Vin and Uint we use Eqs. (116(1 and p7|l to first write 
the site magnetization at zero temperature in terms of 
eigenmodes. Upon taking the ground-state expectation 
value, we have to recall that the vacuum contains zero 
eigenmodes. Hence, 



while 



As a result, 
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(43) 
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-1, iaB. 
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(In order to prove that these are indeed eigenstates, no- 



tice that ^ij — Kii-) These two zero modes do not 

obey the orthogonality relation of Eq. (|28() : they have 
zero hyperbolic norm instead. 



C. Average magnetization per site 

The total staggered magnetization can be written in 
terms of the expectation value of the spin-wave number 
operator: 



^stagg 



\ie^ jeB I 

N 



(39) 
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where we have assumed that the sublattices contain the 
same number of sites: Na = Nb = N/2. As a result, the 



where the sum runs only through eigenmodes with posi- 
tive frequency (the zero modes have been subtracted). 



D. Reduction to a N x N non-Hermitian eigenvalue 
problem 

It is possible to rewrite the 2Nx2N eigenvalue problem 
as two coupled eigenproblems, each one of order N x N 
instead. The non-Hermitian character of the matrices in- 
volved does not change. However, the amount of work for 
numerical computations decreases by a factor of 4 [recall 
that diagonalizing a N x N requires 0{N^) operations]. 

We begin by summing and subtracting Eqs. (??) and 
(??), obtaining 

N 

^^{Kij + Aij){ujn - Vjn) = U!n{Uin + Vin) (46) 



and 
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^{Kij - Aij){Ujn + Vjn) w„(u„ - v^)- (47) 



5 



Multiplying these equations hj K — A and K + A, we 
find the following eigenvalue equations after a simple ma- 
nipulation: 



respectively. Moreover, using these relations and the def- 
inition of 4>''^\ it is straightforward to show that the av- 
erage site magnetization can be written as 



N 



Y}{K~A){K + A)l,{u,^r.- 



(48) 



and 

N 
.7 = 1 



^A){K-A)], 



0- (49) 



Although these equations are in principle decoupled, for 
the purpose of finding the local magnetization they are 
not so, since we are interested in finding mainly v (and 
not u-\-v or u — v alone). We will come back to this point 
later. Equations (|48|l and (|49|l can also be presented in 
more revealing form, namely (here we will drop indices 
to shorten the notation). 



(50) 



where (j)^^^ = u ± v and X — uj^. At this point it is in- 
teresting to notice that the nonzero commutator is the 
cause of non-Hermiticity in the eigenvalue problem. Had 
it been zero, the problem would be become real and sym- 
metric. In fact, it is not difficult to show that 



(51) 



Thus, it is only when all sites have the same number 
or nearest neighbors, i.e., when no dilution is present, 
that the problem becomes real and symmetric (and can 
therefore be solved analytically by a Fourier transform). 
Dilution always makes the eigenproblem non-Hermitian, 
although with real eigenvalues [even if we did not know 
the origin of Eq. (|50|) , it would be easy to prove that all 
eigenvalues A > 0]. 

Let us call M^^^ = if ^ - ± [A, K]. An important 
feature of these matrices is that the left eigenvector (f)^^^ 
of M*^^) is the right eigenvector of M^^^. Thus, if we 
use an algorithm that is capable of finding both the right 
and left eigenvectors of a non-Hermitian matrix, we only 
need to solve the problem for M(+\ for instance. In this 
case, we may say that the 2N x 2N problem has really 
been reduced to N x N. 

In terms of the eigenvectors (t>^^\ the orthogonality 
relations of Eqs. H28() and (|29|l now read 
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and 
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tn Ttjn Tin ^^zm 



= 0, (53) 
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One issue that appears when diagonalizing the problem 
through solving Eq. is that each eigenstate A„ may 
have an eigenvector corresponding to any linear combi- 
nation of the Un and Vn vectors, and not just that u„±u„ 
(that is because each A„ corresponds to at least two eigen- 
frequencies, namely icj, with w„ = Provided that 

there are no other degeneracies, one can sort out which 
combination is generated by noticing the following. Sup- 
pose that 



6*-^^ = c+ (u — v) 



and 



d-) - 



c_(u-fw), (55) 



then, it is easy to see that the normalization conditions 
for both 0*-*' and u,w imply c+c_ = 1. We can then use 
Eqs. and (gU to find that 



and 



{K 
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(56) 



(57) 



These equations provide a way of determining the co- 
efficient c+ (and thus the actual mixing of degenerate 
eigenvectors). For instance, for a given eigenstate. 



1 E. 



N 



/,( + ) 



(K + A) 
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(58) 



Once the c-|_ coefficient has been determined, it is 
straightforward to determine the u„ and Vn vectors cor- 
responding to a given (positive), nondegenerate eigenfre- 
quency w„ in a unique way. If additional degeneracy oc- 
curs, then one needs to introduce more coefficients (and 
consider combinations of all degenerate eigenvectors) in 

Eq. (inni- 



III. NUMERICAL SOLUTION OF THE 
NON-HERMITIAN EIGENPROBLEM 

The numerical solution of the eigenproblem repre- 
sented by Eq. H5U|I requires the full diagonalization of 
at least one real nonsymmetric matrix. However, before 
that, we need to generate the random dilution on a square 
lattice and set the appropriate boundary conditions. An- 
other important point is that we can simplify the diag- 
onalization by breaking the matrix into diagonal blocks, 
each one related to a single disconnected cluster. The di- 
agonalizations can then be carried out separately on each 
block (for each disconnected cluster). Thus, the first task 
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is to reorganize the matrices K and A following a hier- 
archy of disconnected cluster sizes. That involves only 
searching and sorting sites on the lattice (without any 
arithmetic or algebraic manipulation). Moreover, since 
we are interested only in what happens within the largest 
connected cluster (the only one relevant in the thermo- 
dynamic limit and below the percolation threshold)/ we 
can concentrate our numerical effort into the diagonal- 
ization of the matrix block corresponding to that cluster 
alone. 

The first step is to create a square lattice of size N = 
LxL {L being the lateral size of the lattice) with periodic 
boundary conditions in both directions. In order to have 
Na = Nb, we choose L to be an even number. We fix 
the number of holes as the integer part of (1 — p)N and 
randomly distribute them over the lattice with uniform 
probability. 

The second step is to identify all connected clusters 
that exist in the lattice for a given realization of dilu- 
tion. Since most lattices we work with are quite dense 
(relatively few holes), we begin by finding all sites that 
belong to the cluster whose sites are nearest to one of the 
corners of the lattice. Once all sites in that cluster are 
found, they are subtracted from the lattice and the search 
begins again for another cluster. The process stops when 
all sites have been visited and the whole lattice is empty. 
The process of identifying sites for a particular cluster 
is the following. Starting from a fixed site i (the cluster 
seed), we check whether its four neighbors are occupied 
or empty. The occupied ones get the same tag number 
as the first site visited. Then we move on to the next 
side, j -|- 1, and repeat the procedure. We continue until 
we reach the iVth site. 

Along with identifying all sites belonging to each clus- 
ter, we also count then. That allows us to identify imme- 
diately the largest connected cluster in the lattice, whose 
number of sites we call Nc- We set a conversion table 
where the sequential number identifying a site in the 
largest cluster is associated to its coordinate in the orig- 
inal lattice. That allows us to later retrace the compo- 
nents of the eigenvectors of this cluster to their locations 
in the LxL lattice. 

The process of identifying clusters is carried out under 
hard wall boundary conditions. In fact, it is only after the 
largest cluster is found that we force periodic boundary 
conditions. This is the third step. For that purpose, 
we sweep the bottom and top rows, as well as the right 
and left columns of the square lattice, and check whether 
these sites are neighbors of sites belonging to the largest 
cluster once periodic boundary conditions are assumed. 
It turns out that it is easier and faster to do that than 
to search and classify clusters directly from a lattice with 
periodic boundary conditions. 

The fourth step consists of storing the information nec- 
essary to assemble the non-Hermitian matrices of the 
type shown in Eq. (|36|l for the largest cluster only. The 
algorithm is quite fast and allows one to generate and 
find the largest connected cluster in lattices as large as 



L = 100 in less than a second. 

The information generated in the process of identifying 
the largest cluster and its structure is fed into a second 
routine. There, one assembles both the non-Hermitian 
matrices of Eqs. H36|l and (|50|l . We have checked that 
the solution of both the 2Nc x 2Nc and Nc x Nc prob- 
lems provide identical solutions up to several digits for 
a particular realization of the dilution problem at vari- 
ous lattice sizes and dilution fractions. However, only the 
Nc X Nc problem was used to generate the data presented 
here. 

It is important to point out that there exists an alter- 
native formulation of the problem defined by Eq. H13f) . 
using generalized position and momentum operators (see 
Appendix ??). In this formulation one can derive a se- 
quence of transformations that permits the calculation of 
eigenvalues and eigenvectors of the system Hamiltonian 
through the diagonalization of Hermitian matrices alone. 
However, from the computational point of view there is 
no substantial advantage of this approach with respect 
to the non-Hermitian one. 

Since the solution of the non-Hermitian eigenvalue 
problem is less standard than the Hermitian case, we 
provide a description of the method in Appendix ?? 



IV. RESULTS 

We have generated lattices with sizes ranging from L = 
12 to 36 and dilution fractions going from x = to Xc. 
The number of realizations for a given size and dilution 
fraction varied between 500 (nearly clean case) to 1000 
(at the classical percolation threshold). The results of 
the numerical diagonalizations are described below. 

A. Density of states 

The ensemble averaged density of eigenstates as a func- 
tion of frequency for L = 36 is presented in Fig. ?? 
for several dilution fractions and compared to the well- 
known result for the clean case.^^ For a small dilution, 
there is little departure from the clean case, although a 
small structure is already visible at around uj/JS = 3. 
As the dilution increases, a peak and an edge develop at 
around this frequency. Notice that the overall trend is 
a decrease in the number of high-frequency modes, with 
the proportional increase in the number of low-frequency 
ones. Close to the percolation threshold, another struc- 
ture appears at around to/ JS — 2. Thus, we see that the 
effect of dilution is to shift spectral weight from high to 
low frequencies in a nonuniform way. This tendency had 
also been observed in Rcf. 9 for small dilution. 

Two additional very sharp peaks (not shown in Fig. 
??) also exist in the density of states as the dilution in- 
creases. They occur at frequencies uj/JS =1,2 and cor- 
respond to configurations where Vin = for all sites in 
the cluster, while Um = for all but two sites. Their 
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FIG. 1: Average density of states as a function of energy 
for four different site dilution fractions: x — (solid line), 
X = 0.1 Xc (dotted line), x = 0.5 (dotted-dashed line), and 
X = Xc (dashed line). Notice the two structures just below 
uj/JS = 2 and 3. The latter is already visible at a; = 0.1 a;^, 
while the former becomes prominent only for x > QAxc- The 
data was obtained from 500 to 1000 realizations of a 36 x 36 
lattice. 



-0- 



(a) 




O- 
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FIG. 2: The two dangling structures that occur in the lo/ J = 
1 (a) and 2 (b) eigenstates. The filled circles indicate iti„ ^ 0, 
while the empty circles have Uin = 0. All sites in these states 
have Vin = 0, thus they do not contribute to the staggered 
magnetization. 



typical spatial structures are shown in Fig. ??. Since 
Vin = for all sites, these states do not contribute to the 
quantum corrections to the staggered magnetization. 

For the clean case, it is simple to verify (based on the 
exact diagonalization of the problem) that the low fre- 
quency modes provide the largest contribution to Smz- 
For the dilute lattice, the same is true, as can be seen in 
Fig. ??, where we have plotted 
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FIG. 3; Staggered magnetization per unit of magnetic site as 
a function of energy, 5mz(tj), under the same conditions of 
Fig. ??. 



B. Inverse participation ratio 

The nature of the eigenstates also changes as the dilu- 
tion increases. The best way to characterize the nature 
of the states is through the return probability, that is, 
the probability that after some very long time a particle, 
moving in the percolating lattice, will return to its origi- 
nating point. The return probability can be expressed in 
terms of the inverse participation ratio (IPR).^^ Here, we 
use a definition of the IPR involving the eigenvector com- 
ponent related to the quantum fluctuation corrections to 
the magnetization, namely. 



T,n=2 H'^ 



where 



T,n=2 ^i"^ - 



2^1=1 ' 



(60) 



(61) 



In Fig. ?? we show the IPR as a function of energy for 
three lattice sizes. According to its definition, the IPR 
for extended states decrease as the system size increases, 
while for localized states the IPR is insensitive to any 
size variation. These trends are clearly visible in Fig. 
??, namely, states are mostly extended when dilution is 
small and tend to localize as one gets closer to the perco- 
lation threshold. For intermediate dilution [Fig. ??(c)], 
we see that the states close to uj/JS = 3 are strongly 
localized while the remaining states are quite extended. 
As expected, the low frequency states tend to remain ex- 
tended up to strong dilution. 



As a result, we see that the transference of eigenmodes 
from high to low frequencies is the mechanism by which 
quantum fluctuations are enhanced as the dilution in- 
creases. 



C. Average magnetization 

In the thermodynamic limit, the magnitude of the stag- 
gered magnetization per unit of lattice site can be 
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CO/JS CO/JS <B/JS 



FIG. 4: Average inverse participation ratio, as defined in Eq. 
1601 1. as a function of energy for different lattice sizes and 
dilution fractions: (a) x = 0.1 Xc, (b) x = 0.5 ajc, and (c) x = 
Xc- As the dilution increases, states become more localized, 
beginning with those located in the high-energy part of the 
spectrum. Each curve shown corresponds to an average over 
500 to 1000 realizations. 



written as 




{S-5m,), (62) 



where Nm is the total number of occupied sites in the lat- 
tice {Nm = pN). While the first factor on the right-hand 
side of Eq. (|62|l is purely classical, the second factor is 
purely quantum, namely, it measures how quantum fluc- 
tuations reduce magnetic order. Thus, since Nc vanishes 
at Xc, a quantum critical point below Xc can only exist 
if, for some x* < Xc, we find Sruz = S. If, on the con- 
trary, SiTiz < S a,t X = Xc, then the order is only lost 
at the percolation threshold and the transition is essen- 
tially classical. It is important to have in mind that the 
linear spin- wave approximation is well defined only when 
6mf < S, for all i = 1, . . . , iV. When Sml « S" at a 
large number of sites, the approximation is not necessar- 
ily quantitatively correct. 

Figure ?? shows the average staggered magnetization 
per unit of magnetic site as a function of dilution fraction, 
mz{x), when S = 1/2. The points were obtained after 
finite-size scaling the ensemble averaged data taken from 
12 different lattices sizes. As a consistency check, we have 
also calculated the staggered magnetization for the clean 
case {x = 0, no ensemble average) with the same numer- 
ical procedure. We have found that mz{0) ~ 0.303, con- 
sistent with values obtained by other methods. Thus, 
at least at low dilution, the spin-wave approximation is 
quite accurate. 

For comparison, we also show experimental data ob- 
tained for LCMO from both neutron scattering^ and nu- 
clear quadrupole resonance (NQR)f^ as well as the re- 
sult of QMC simulations of the dilute Heisenberg AFM 
in a square lattice>^ One can see that our simulations, 
based on the spin-wave approximation, capture the main 



FIG. 5: Average staggered magnetization per unit of mag- 
netic site. Results of the S — ^ SWT simulations after finite- 
size scaling (circles) are compared with neutron scattering 
data (squares)^ NQR data (diamonds) r" and the fit to the 
QMC data from Ref. @ (solid line). Also shown is the occupa- 
tion fraction of the largest connected cluster, Nc/Nm, (dashed 
line) and the analytical result from the calculation of Ref. |^ 
(dashed-dotted line). Inset: the quantum fluctuation contri- 
bution to the staggered magnetization for different values of 
S. 



features of the experimental data, namely, a progressive 
decrease of the staggered magnetization up to the clas- 
sical percolation threshold. At a dilution fraction very 
close to Xc, our simulations indicate that the staggered 
magnetization should vanish. The inset in Fig. ?? shows 
that the vanishing of the staggered magnetization occurs 
because Sm^ goes to 1/2 very close to the classical transi- 
tion point. Thus, the same effect would not arise had we 
used S > 1/2. The QMC simulations, on the other hand, 
predicts dniz < 1/2 at x = Xc, thus indicating that the 
transition is purely classical. The relatively small num- 
ber of experimental points and the large error bars near 
the percolation threshold do not allow for an adequate 
distinction between a classical and a quantum transition 
for LCMO. 

The discrepancy between our result and the QMC 
simulations for the staggered magnetization close to Xc 
should be seen as an indication that, while qualitatively 
correct, our approach fails quantitatively when the or- 
der parameter magnitude is significantly reduced locally. 
This is expected if we recall the assumption used in the 
derivation of Eq. H10|l . Nevertheless, the spin- wave ap- 
proximation, having access to low-lying excited states 
and wave functions, allows us to understand in more de- 
tail, at least qualitatively, how the suppression of order 
due to quantum fluctuations takes place upon dilution. 
This is not the case for the QMC simulations. In fact it 
is surprising that our calculations seem to agree with the 
experimental data better than the QMC. This can be un- 
derstood by the fact that the experimental system may 
contain extra oxygen atoms that introduce holes in the 
Cu02 planes, as well as next-nearest neighboring inter- 
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FIG. 6: The density plot of local demagnetization Smf for 
a 32 X 32 lattice with periodic boundary conditions at a:: ~ 
Xc- Only sites belonging to the largest connected cluster are 
shown. The weak links present in the cluster are indicated by 
bullets. 



actions, that frustrate the AFM state and also introduce 
larger quantum fluctuations. The latter are captured by 
the overestimation produced in the linear spin-wave the- 
ory. In fact, it is known that La2Cu04, has a nonzero 
frustrating next-nearest neighbor coupling.^ That effec- 
tively decreases the spin per site to a value smaller than 
1 /2, possibly bringing LCMO closer to a quantum critical 
point than the pure 5 = 1/2 Heisenberg AFM. 



D. Local fluctuations 

We now turn to the question of local fluctuations. We 
have so far discussed the site-averaged demagnetization 
Sniz and used the criterion that it must be smaller than 
S for the order to persist. However, one could argue that 
some sites may have particularly large fluctuations; if 
these large fluctuations take place exactly at weak links of 
the largest connected cluster backbone, then they could 
be responsible for earlier destruction of the long-range 
order. We have numerical evidence that this is unlikely, 
although the relatively small size of our lattices does not 
allow us to be conclusive. In Fig. ?? we show an intensity 
plot of the local quantum fluctuations Srui in the largest 
connected cluster, very close to the percolation threshold, 
for a typical realization of a L = 32. Notice that the 
largest fluctuations tend to appear only along deadends 
or dangling structures, and not in the links connecting 
blobs of the cluster backbone. The same trend is seen in 
all realizations that we have inspected. 



E. Excited states 

Equation (|13(l represents a system of N coupled har- 
monic oscillators with A = K + A and B = K — A 
being the spring constant and the inverse mass tensors, 
respectively (for an alternative description of the prob- 
lem in terms of position and momentum operators, see 
Appendix ??). For the simple model of elastic vibrations 
in a lattice, when qi — [ai a\ )/V2 has the meaning 
of a displacement of the ith atom around its equilibrium 
position, it is well known that the Hamiltonian of such 
a system can be mapped onto a problem of diffusion in 
a disordered lattice.^ Analytical results for the related 
diffusive problem, as well as numerical simulations with 
large systems, allow one to understand several properties 
of the diluted vibrational model, such as the density of 
states (DOS) and the dynamical structure factor. Per- 
haps one of the most distinctive features is the existence 
of strongly localized excitations named fractonsi^S For 
dilution fractions x < Xc, these excitations appear in the 
high-frequency portion of the spectrum, beyond a cer- 
tain crossover scale ujc, while the low- frequency part is 
dominated by acoustic phononlike, nearly extended, ex- 
citations. At exactly x = Xc^ the systems becomes a 
fractal and the fractons take over the entire spectrum^ 
Scaling considerations, as well as numerical simulations, 
have shown that, for a square lattice, the DOS behaves 
as 



(63) 



The crossover frequency depends critically on the dilu- 
tion: UJc (X {xc - x)^, where D « 91/48. This IS an 
important result because, for the elastic vibration prob- 
lem, it is also possible to show that the quantum fluc- 
tuation corrections to the classical order parameter obey 
the relation 



dm 



1 

N 



N 

E 



1 j V^o o 



duj 



(64) 



where {w„} are the nonzero eigenvalues of the corre- 
sponding Hamiltonian. Therefore, Sm remains finite be- 
low and at the percolation threshold, indicating that 
quantum fluctuations are likely not sufficiently enhanced 
by the dilution to destroy the existent long-range orderii^i 
We will get back to this point in Sec. ??. 

In order to investigate the existence of such fractons in 
the dilute 2D Heisenberg AFM, as well as to clarify the 
nature of its low-lying excitations, we have calculated the 
the dynamical structure factor, 

^ gjq-(R,-Rj) 

x{S+{0)S^{t) + S-{0)S+{t)). (65) 

By using Eqs. (jSJ, ®, and (fH^ . we can express 5(q, lj) 
in terms of Fourier transformations of the site-dependent 



iS(q, w) = dte 
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Bogoliubov coefRcients Um and Vin- We find 



2S 



5{iO-UJn) Un{^) "n(-q) 



-q) 



^n(q) ^rf(-q) 



+f)f(q) ^;f(-q)+5^(q) S^(-q) 



+u^(q) w^?(-q) + 



(66) 



where the partial terms involving u^'^ are given by the 
Fourier transformation of u„, namely, 



(q) 



(67) 



and analogously for w^''^. Notice that the sum over sites 
runs only over one of the sublattices, A or depending 
on the particular term. Thus, only four two-dimensional 
Fourier transformation are required in order to evaluate 
5(q,w). 

We have computed numerically the Fourier transfor- 
mations and calculated the average dynamical structure 
factor for lattices of size L = 32 at several dilution frac- 
tions. Only sites within the largest connected cluster 
were taken into account. Averages were performed over 
50 realizations for each case. The results are presented in 
the form of intensity plots in Fig. ??. To provide a bet- 
ter contrast, we have rescaled (iS(q, a;)) by the function, 
f{uj) — X]q('5(q; ^))- Only the data along two partic- 
ular directions in momentum space are shown, namely 
Q = Qx = Qy and q = Qx, Qy — 0, with < g < tt (the lat- 
tice spacing is taken to be unit). For small dilution, the 
structure factor resembles closely that of the clean case, 
with some small broadening of the magnon branch due 
to the weak destruction of translation invariance. How- 
ever, particularly along the qy = direction, one can 
already notice a small hump at around uj/JS = 3, con- 
sistent with the peak-and-edge structure seen in Fig. ??. 
This feature becomes more prominent with increasing di- 
lution. For dilution fractions larger than O.Gxc, another 
hump becomes visible at around lu/JS = 2, again con- 
sistent with the feature observed in Fig. ?? at the same 
frequency. Close to the percolation threshold, there ex- 
ist three clear broad branches in the spectrum. While 
the dispersion of the high-frequency branch at uj/JS > 3 
is hardly affected by the dilution, the opposite occurs 
with the low-frequency one, at lu/JS < 2, where the 
slope (magnon velocity) decreases with increasing dilu- 
tion. We interpret the progressing breaking and bending 
of the magnon branch as the system becomes more di- 
luted to the appearance of fractons. At x — Xc the exci- 
tation spectrum has little resemblance to that of x = 
and even the long wavelength part is strongly modified. 
In between these two limits, there is a crossover from a 
magnon-dominated to a fracton-dominated spectrum. A 
three-branch structure for the spectral function in the 
spin- wave approximation was also found in Ref. 9. How- 
ever, the positioning of those branches and their relative 
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FIG. 7: Gray plots of the rescaled average dynamical struc- 
ture factor sliced along two distinct directions in momentum 
space: qx = qy = q [(a), (c), and (e)] and qx = q,qy =G [(b), 
(d), and (f)]. Plots at different dilution fractions do not nec- 
essarily have the same arbitrary scale for the gray intensity. 
The most prominent branches are marked with dashed lines. 
Averaging for each plot was performed over 50 realizations of 
dilution. 



intensity were different from what we observed in our nu- 
merical solution. We believe the cause of this discrepancy 
is the limited range of applicability of the perturbative 
treatment, which is expected to be accurate only in the 
weak dilution regime. 

The gradual appearance, broadening, and motion of 
the branches are better represented in Fig. ??, where the 
rescaled average dynamical structure factor is shown as a 
function of frequency for a fixed momentum qx = OA/n 
and qy — 0. The two-peak structure observed in our 
numerical data has some resemblance to the results of 
inelastic neutron scattering performed by Uemura and 
Birgeneau for the compound Mna;Zni_a;F2fi2. whose mag- 
netic properties are described by the three-dimensional 
site-diluted random Heisenberg model with S = 5/2. 
These authors observed two broad peaks which were asso- 
ciated with magnons (low-energy, extended) and fractons 
(high-energy, localized) excitations. The relative ampli- 
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FIG. 8: Slice of the rescaled average dynamical structure fac- 
tor for Qx = OA/pi and = for L = 32 and different 
dilution fractions. 



as the high-frequency states become more localized with 
increasing dilution. This feature is masked in Fig. ?? 
by the frequency-dependent rescaling, but is clear from 
Fig. ??. The high-frequency, high momentum modes 
are much less dispersive than the low-frequency, low mo- 
mentum ones. In fact, particularly below lu/ JS ~ 3, the 
magnon branch is q independent (i.e., broad in momen- 
tum) and thus strongly localized. 

The low-frequency modes are likely not fully ballis- 
tic (coherent), given the large broadening seen in the 
structure factor, but should rather have a diffusive prop- 
agation at large scales. As to whether they remain two 
dimensional or gain a lower-dimensional character, our 
data is not conclusive. 



V. UPPER BOUND FOR QUANTUM 
FLUCTUATIONS IN BOSONIC SYSTEMS 



tude, width, and dispersion relations of these excitations 
were measured and were found consistent with the theo- 
retical predictions by Orbach and Yu^^ and Aharony and 
co-workers i2S The scaling theory of the latter predicts 
that a two-peak spectrum appears at momenta q ~ 
where ^ is the percolation correlation length (i.e., the typ- 
ical linear size of the disconnected clusters when x < Xc)- 
The relatively small size of the L = 32 lattice does not 
allow us to make a similar quantitative comparison be- 
tween these theories and our numerical results. 

While we do observe a two-peak structure, we also 
see a third peak at sufficiently large diluation, although 
quite weak. This additional peak may be characteristic of 
square lattices; in the body-centered-cubic Mna;Zni_a;F2, 
there might exist a different multipeak structure. More- 
over, the large broadening and limited energy resolution 
of the neutron-scattering experiments may have made 
such features unobservable in the data of Ref. 0. How- 
ever, one aspect which seems different between the exper- 
imental data and our numerical simulations is the way 
the spectral weight is transfered between magnons and 
fractons as a function of dilution. While Uemura and 
Birgeneau see an increase (decrease) of high-frequency 
fractons (low-frequency magnons), we see the opposite: 
The low-frequency portion of the spectrum becomes more 
intense. 

The averaging procedure recovers, to some extent, the 
translation invariance broken by the dilution. This al- 
lowed us to identify the large momentum part of the 
branches which, otherwise, would not be visible. How- 
ever, the low momentum part is clearly visible even when 
no averaging is performed on the data (not shown) . This 
fact, together with the strong dispersion of the lower 
energy branches, suggests that long wavelength propa- 
gating magnons are present for dilution fractions below 
the percolation threshold, and lose out to fractons at 
X = Xc- At T = 0, these modes are responsible for 
the quantum fluctuations that bring down the long range 
order. Their spectral weight becomes more important 



The possibility of a classification of quantum random 
systems into universality classes is an important theo- 
retical problem with relevance to experiments. Random 
matrix classification schemes introduced by Wigner lim- 
ited the classes of problems initially to the orthogonal, 
unitary, and symplectic classes^^ but these were later 
extended to encompass the classes of chiral models^'* and 
models where the particle number is not a conserved 
quantity.^^ While most of the random matrix problems 
are related to fcrmionic spectra, there is renewed interest 
in the problem of bosonic random matrix theoryj2£ Of 
particular importance is its application to the problem 
of diluted quantum magnets since these systems, in the 
limit of large spin S, can be approximately described, via 
linear spin-wave theory, as a problem of noninteracting 
bosons 

One of the main differences between fermions and 
bosons is that, in addition to the symmetries of the un- 
derlying Hamiltonian, one must ensure that the bosonic 
spectrum is semipositive definite; this stability condition 
is not an issue in fermionic systems. However, it is au- 
tomatically satisfied when the disorder is caused by site 
dilution. 

Let us concentrate our discussion on systems whose 
Hamiltonian can be mapped onto a set of N coupled 
harmonic oscillators of the following kind: 

1 ^ 

= 2 II 1i + Pj ) • (68) 

Here, N is the total number of sites of the square lattice, 
while qi and pi represent generalized position and mo- 
mentum operators at a lattice site i, respectively, such 
that [qi,Pj] = iSij, with i,j — 1,...,N. The N x N 
matrices A and B are real, symmetric, and semiposi- 
tive; in the most general case, they do not commute. 
The magnitude of quantum fluctuations are character- 
ized by the mean-square deviation of the average value of 
these operators in the ground state: q^ = J^iilf)/-^ ''^^'^ 
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= ^i{Pi)/N. As mentioned in Sec. ??, if quantum 
fluctuations are unbounded at the percolating regime, 
then these mean values diverge and LRO is not possi- 
ble, implying that order has to be destroyed before Xc 
is reached; that is, a quantum critical point should exist 
before the percolation threshold. Here we show that this 
is not the case and that LRO can persist up to and in- 
cluding Xc- As for the spin Hamiltonians approximated 
by Eq. the exact disappearance of magnetic LRO 

may also depend on the spin value. 

The eigenfrequencies of the harmonic oscillators of Eq. 
(|68|l can be shown to be those of 



ni = AB or nj^ = BA. 



(69) 



Since, in general, [A, B] is nonzero, the matrices and 
are non-Hermitian and distinct. However, it is simple 
to show that they have exactly the same real eigenvalues. 



for n — 1, 



,N. 



This is also true for the matri- 
^i/2^^i/2_ Therefore, 
the fluctuations of the position and momentum, averaged 
over all sites, can be written as^l 



ces ill = B^/^AB^I"^ and ill 



f = J_tr (Bilc-^) < - At, 



P 



— tr(Anc ')<—i^, 
2N \ ^ y - 2 ' 



(70) 
(71) 



where and 6, denote the maximum eigenvalues of the 
semipositive definite matrices A and B, respectively, and 
we used that 



K = — tr f2^^ = — tr Clr,^ = 



N 



N 



duj p{lu)uj-^, (72) 



since He and Clc share the same eigenvalues (wmax being 
the largest), with spectral density p{uj)- Thus, the finite- 
ness of the quantum fluctuations reduces to the problem 
of the convergence of the integral in Eq. (??) (quantum 
mechanics requires that (f > 1/4). Notice that the 
quantum fluctuation correction to the order parameter 
per imit of lattice site can be written as 



(73) 



where a, (3, and 7 are constants that depend on the par- 
ticular model and order parameter under consideration. 

We now turn to applying these general results to spe- 
cific bosonic models based on Heisenberg Hamiltonians. 



A. 0(2) model 

This model is realized, for example, in an array of 
Josephson junctions, where the variables qi correspond 
to the linearized phase of the superconductor order 
parameter^ '^'^^ The matrices in Eq. H68() take the sim- 
ple forms Aij — Kij — Ay and Bij = {U/J)Sij, where 
Aij = 1(0) when the nearest-neighboring sites i,j are 
both occupied (otherwise) and Kij = Sij "^iLi^ih 



Kii counts the number of nearest neighbors of site i. 
Here, J denotes the Josephson coupling and U is the is- 
land charging energy (usually, J ^ U). The same struc- 
ture occurs in the case of vibrations in a diluted latticed 
For the 0(2) model, the frequency eigenvalues are ob- 
tained from il^ = (U/J) A. The problem of determining 
the density of states of the connectivity matrix A can 
be mapped onto a problem of diffusion in a disordered 
latticei^ It can be shown that the density of states of 
il goes p{uj) (X oj'^ Up to the percolation threshold, 
d* — d — 2, while d* = J « 4/3 exactly at x = a;c. There- 
fore, the site-averaged fluctuations are bounded, and the 
linear approximation (as well as order) is maintained as 
long the ratio J/U is large enough. 



B. Heisenberg antiferromagnetic model 



When we take Jx — Jy — Jz, the matrices take the 



form Aij = Kij - 



- Aij and Bij — Kij 



- Aij . It is useful to 



define a matrix A = A that has the effect of changing 
the signs of qi and pi for all sites i in one of the two 
sublattices. Thus, B = AAA, and = AAAA. Notice 
that r = A^, is not semipositive definite (in fact it is 
non-Hermitian), but has real eigenvalues with the same 
magnitude as those of il. 

In the 0(2) model, because Sl^ cx A, one can directly 
relate the energy eigenvalues of il to those of the matrix 
A. In the AFM case, we need to obtain the density of 
states of r = AA, which is not simply related to that of 
A (since [A, ^] ^ in general for diluted systems). This 
nontrivial relation between the eigenvalues of il and A 
is generally present in bosonic problems. For example, 
a similar problem also appears in the work of Gurarie 
and Chalker in the relation between their stiffness and 
frequency eigenvalues.^^ 

The problem of determining the density of states of 
r for a random dilution problem is one of the interest- 
ing open questions related to the important difference 
between random fermionic and bosonic systems. In a 
fermionic problem this question would be already an- 
swered by matching the symmetries of F to the Cartan 
classification table. However, one cannot substitute the 
matrix A by an arbitrary random matrix with similar 
symmetries, since that would violate the semi-positive 
definite constraint. In this work we do not attempt to 
analytically resolve this problem; however, we find nu- 
merical evidence that the density of states of F^ follows 
that of A at low energies. 

Our simulations show that the site-averaged fluctua- 
tions Smf are bounded, both below and at the percola- 
tion threshold. This means that order should exist up to 
and including the critical dilution Xc, as long as Srui is 
small compared to the value of the spin S. Recalling Eq. 
(|62|l . we conclude that there could be a minimum value 
Smin for the spin, below which there is a quantum phase 
transition for dilutions x < Xc- An effective local spin 
smaller than 1/2 can be realized in a bilayer system with 
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antiferromagnetic interlayer couplingi^^ 



C. XXZ model 

In this case, Jx — Jy ^ Jz- We then have Aij — 
Kij — and Bij = Kij — 7 , where 7 measures the 
anisotropy. Alternatively, we can write B = (l + 7)/2A + 
(l-7)/2AAA, and so Vl^ = (l + 7)/2A2 + (l-7)/2AAAA 
(notice that for 7 ^ —1 we have the same problem as 
the AFM). The analysis is similar to the previous two 
cases. The amount of fluctuations is controlled by the 
anisotropy, and it can be shown to be bounded if the 
density of states follows that of the AFM case. 

VI. DISCUSSION AND CONCLUSIONS 

In this work we have studied the role played by site di- 
lution in enhancing quantum fluctuations in the ground 
state of the Heisenberg antiferromagnet in a square lat- 
tice. Using the linear spin-wave approximation for this 
model, we have performed exact numerical diagonaliza- 
tions for lattices up to 36 x 36, with dilution fractions 
going from the clean limit to the classical percolation 
threshold. Our results indicate a progressive, nonuniform 
shift of spectral weight in the spin-wave excitation spec- 
trum from high and to low frequency as the dilution in- 
creases, with the high-frequency part becoming more lo- 
calized. The higher density of low-frequency, long wave- 
length excitations leads to strong quantum fluctuations 
and a decrease in the magnitude of the staggered magne- 
tization. For dilutions very close to the classical percola- 
tion threshold, we have found that quantum fluctuations 
are sufficiently strong to nearly match the clean-limit 
magnitude of the magnetization when S — 1/2, but not 
for higher spins. This is consistent with recent neutron- 
scattering experiments with the S = 1/2 dilute Heisen- 
berg antiferromagnet La2Cui_a;(Zn,Mg)a;04, which show 
that long-range order disappears at around the classical 
percolation threshold. However, quantum Monte Carlo 
simulations suggest that quantum fluctuations should re- 
main small and that the destruction of long-range order 
is controlled only by the disappearance of the infinite 
connected cluster. We understand this discrepancy be- 
tween the quantum Monte Carlo results and the linear 
spin-wave theory near the classical percolation threshold 
as an indication that the latter has its validity limited, 
as the magnitude of the order parameter is too small. 

While perhaps not quantitative-accurate, our simula- 
tions do allow us to probe into the nature of the low- 
lying excited states of the dilute antiferromagnet. We 
observe two clear humps in the density of states at fre- 
quencies (jj/JS = 2 and 3. By calculating the ensemble- 
averaged dynamical structure factor, we were able to as- 
sociate the appearance of these humps with the breaking 
of the clean-limit magnon branch into three distinct but 
broad branches. The new branches tend to be strongly 



localized (nondispersive) at high frequencies and have a 
diffusive rather than ballistic nature at low frequencies. 
In the literature, the multipeak structure had been asso- 
ciated with the appearance of fractons in the excitation 
spectrum as the dilution increases. From our simula- 
tions, it seems that the picture is somewhat more com- 
plex. Besides the overall broadening, the position of the 
high-frequency branch remains close to the corresponding 
portion of original clean-limit magnon branch, while the 
magnon velocity (the dio/dq slope) in the low- frequency 
branch is continuously reduced with increasing dilution. 
Therefore, it appears that the fracton character also con- 
taminates the low-frequency branch. However, the lack 
of resolution due to the finite size of our lattices does 
not allow for a conclusive picture. We did not attempt, 
however, to study fracton states which can possibly occur 
above the percolation thresholdiSSiSS 

It is important to remark that, in principle, due to An- 
derson localization in two dimensions, we expect that in 
the infinite system all excitations should in fact be local- 
ized for any finite dilution. In order to probe more care- 
fully strong localization and the consequent exponential 
decay in real space, we need not only much larger lat- 
tices, but also to calculate two-point correlators, which 
goes beyond the applicability of linear spin-wave approx- 
imation. For that same reason, we were not able to eval- 
uate quantities such as the spin stiffness, which involves 
matrix elements of higher than bilinear operators. 

We also studied the question of local quantum fluctua- 
tions as a way to destroy long-range order. For finite-size 
lattices we found that the weak links do not show strong 
quantum renormalizations. That provides some indica- 
tion that local quantum fluctuations may not be sufficient 
to change the dominance of the classical percolation pic- 
ture. 

Using a more general analytical formulation, we have 
argued that there exists an upper bound for the quan- 
tum fluctuations in any model with a classically ordered 
ground state whose Hamiltonians can be mapped onto 
that of a system of coupled harmonic oscillators. The 
amount of quantum fluctuations depends directly on the 
low-energy behavior of the density of states of the asso- 
ciated bosonic problems. Our exact diagonalization of 
the linear spin- wave Hamiltonian on a percolating lattice 
led us to identify the value of the upper bound for one 
particular type of model and can readily be used to find 
similar values for any other bosonic model. This could 
be used to study a large class of spin Hamiltonians that 
can be bosonized in the ordered phase. 
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APPENDIX A: FORMULATION IN TERMS OF 
COUPLED HARMONIC OSCILLATOR 



alternatively, to set q'g = 0): {q'}N W}n-i and 
{p'}n {p'}n-i- Also, [AVxAf [A']n-ixN-i and 
[b]jvxAf [b]7V-ixAf-i- Notice that now all bk > 0, 
k — 2, ... TV. Thus, in this new space, we can perform 
the following rescaling: 



" = b-i/2g' and p" = h^/^p', 



(A8) 



The spin-wave Hamiltonian of Eq. H13() can be rep- 
resented in terms of position and momentum operators. 
In this language, it becomes more transparent that the 
problem of finding the eigenvalues and eigenvectors of 
the Hamiltonian can be solved by the diagonalization of 
two real symmetric matrices, an alternative to the non- 
Hermitian eigenvalue formulation of Sec. ??. 

Let us perform the following operator transformation: 



ai + a. 



and Pi = 



iV2 



(Al) 



for all i = 1,. ..,N. Notice that the operators qi and 
Pi obey the canonical position-momentum commutation 
relations. It is convenient to adopt a matrix formulation 
for the problem, namely. 



H^:l^{q^Aq+p^Bp) 



(A2) 



where x = {xi}i^i,,,N and p = {pi}i=i...N denote vec- 
tors of position and momentum operators, while [AJ^ = 
Kij + Aij, and [B].y — Kij — Ay. The quantum fluctu- 
ation correction to the sublattice magnetization can be 
written as 

Sm=^tr[{qq'^) + {pp^)]~^. (A3) 

We can diagonalize B through an orthogonal transfor- 
mation U, 

U^B U = b (diagonal) (A4) 
and define new coordinates such that 

q = Uq and p = V p' . (A5) 
Defining A' = IJ-^A U, we then have that 

H^:^{q'^A'q'+p'^hp') (A6) 

and 

^™-^[tr((<z'9'"))+tr((p'y-))]-i. (A7) 

It is not difficult to prove that all elements in the di- 
agonal of b are positive except one, bi, which is zero. 
In order to eliminate this zero mode, we subtract the 
corresponding row or line in all vectors and matrices, 
which amounts to a reduction in the Hilbert space (or. 



which allows us to write 
JS 



H=^{q"^A"q"+p''^p"), (A9) 



with A" = bi/2A' bi/2 as well as 



1 



Sm = ^ [tr(b +tr(b-l {p" p"^)) 



Hp?)]--,, 



(AlO) 



where the last term within the square brackets reflects 
the existence of a zero mode. It is useful now to return 
to the site basis by carrying out the inverse rotation. 



such that 



q'" = Vq" and p'" = Up", 



H=^{q"'^Cq"'+p'"^p"'] 



(All) 



(A12) 



where C = U A" U"^ is the new connectivity matrix. If 
we define 

B^/2^Ubi/2u^ and B-^/^ = U b-^/^ (A13) 

within the subspace — 1 x — 1 where no zero mode 
is present, we can write the connectivity matrix as 



C = B1/2ab1/2. 



(A14) 



Also, notice that the inverse rotation does not change the 
expression of the sublattice magnetization, 

Sm ^ ^ [tr(b-i +tr(b (p'>"'^)) 



(A15) 



We can now perform the last operation, namely, the 
diagonalization of the connectivity matrix. 



= Vq'" and = Vp'" 



yielding 



where 



V^CV^c (diagonal). 



(A16) 



(A17) 



(A18) 
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with all c„ > 0, n = 1, . . . — 1. We finally arrive to a 
system of decoupled harmonic oscillators. The quantum 
fluctuation part of the magnetization becomes 



+tr(bV^ (P™P'"^)V) + H2) 
However, we know that 



and 



-(A19) 



(A20) 



(A21) 



Therefore, 



Sm = 



1 

iiV 



tr(b-iV^c-i/2v) +tr(bV^ci/2V 
1 

2' 



(A22) 



It is interesting to notice that, since A = OBO, we 
have that 



where 



n = bi/2qbi/2 



(A23) 



(A24) 



Thus, it is easy to see that if C were a positive matrix, 
then we would be able to write c^/^ V = C^/^ ~ fl. 
That would allow us to simplify the expression for the 
sublattice magnetization a step further. However, the 
connectivity matrix is not necessarily positive. 



APPENDIX B: NUMERICAL 
DIAGONALIZATION OF THE NON-HERMITIAN 
MATRIX 



The diagonalization of the non-Hermitian matrices 
consists of five steps. First, the real (but asymmet- 
ric) matrix A/^^^ is reduced to an upper Hessenberg 
form through an orthogonal transformation, namely, A = 
QM^+^Q^. This is done by using the LAPACK subrou- 
tines DGEHRD and DORGHR. Second, we use the LA- 
PACK subroutine DHSEQR to perform the Schur fac- 
torization of the Hessenberg matrix: A = ZTZ^ . That 
allows us to obtain the eigenvalues and the Schur vectors, 
which are contained in the orthogonal matrix Z . Third, 
using another LAPACK subroutine, DTREVC, we ex- 
tract from Z both right and left eigenvectors of M(+^ 
In the fourth step we renormalize all eigenvectors {(/>'-^^} 
such that they satisfy Eq. (|52|) [the condition in Eq. 
(|53|l is automatically satisfied], sort the eigenvalues {A„} 
in ascending order, and extract the zero mode from the 
spectrum. For some realizations, the lowest eigenvalues 
next to the zero mode cannot be distinguished from the 
zero mode itself and are therefore neglected. Finally, the 
correct linear combination of (^^^^ that provides the cor- 
rect u„ and Vn for each positive eigenfrequency = \f\i 
is obtained according to the algorithm presented in Sec. 
?? 
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